<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of v_psycest</title>
  <meta name="keywords" content="v_psycest">
  <meta name="description" content="V_PSYCEST estimate multiple psychometric functions">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../index.html">Home</a> &gt;  <a href="index.html">v_mfiles</a> &gt; v_psycest.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for v_mfiles&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->

<h1>v_psycest
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>V_PSYCEST estimate multiple psychometric functions</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>function [xx,ii,m,v,mr,vr]=v_psycest(iq,x,r,xp,lf) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre class="comment">V_PSYCEST estimate multiple psychometric functions

 Usage: [xx,ii,m,v]=v_psycest(-n,p,q,xp,lf) % initialize n models
        [xx,ii,m,v]=v_psycest(i,x,r)     % supply a trial result to v_psycest
        [xx,ii,m,v]=v_psycest(i,x,r,o)   % supply a trial result to v_psycest and plot
                    v_psycest(i,o)       % plot pdf of model i with plot options o
        [xx,ii,m,v,mr,vr]=v_psycest(i,o) % Determine robust outputs in addition expcluding outliers [takes longer]
              [p,q,msr]=psychest(0)    % output model parameters (or print them if no outputs)

 Inputs (see usage examples for argument positions):
         -n        minus the number of models
          p(:,n)   parameters for each model
                      1  thresh [0.75]
                      2  miss [0.04]
                      3  guess [0.1]
                      4  SNR min [-20]
                      5  SNR max [20]
                      6  Slope min (but over-ridden if &lt;q.sl) [0]
                      7  slope max [0.5]
          q(:)     parameters common to all models (vector or struct)
                      1  q.nx  number of SNR values in pdf [40]
                      2  q.ns  number of slope values in pdf [21]
                      3  q.nh  number of probe SNR values to evaluate [30]
                      4  q.cs  weighting of slope relative to SRT in cost function [1]
                      5  q.dh  minimum step size in dB for probe SNRs [0.2]
                      6  q.sl  min slope at threshold (must be &gt;0) [0.005]
                      7  q.kp  number of std deviations of the pdfs to keep [4]
                      8  q.hg  amount to grow expected gains in ni trials [1.3]
                      9  q.cf  cost function: 1=variance, 2=v_entropy [2]
                     10  q.pm  psychometric model: 1=logistic, 2=cumulative gaussian [1]
                     11  q.lg  use log slope in pdf: 0=no, 1=yes [1]
                     12  q.pp  Number of prior standard deviations in initial semi-range [1]
                     13  q.pf  Probability floor (integrated over entire grid) [0.0001]
                     14  q.ts  Number of std devs to explore [2]
                     15  q.dp  Maximum probe SNR shift (dB) [10]
                     16  q.it  Grid interpolation threshold [0.5]
                     17  q.at  Axis change threshold [0.1]
                     18  q.la  Look 2-ahead when choosing probe [1]
                     19  q.op  Outlier probability [0.01]
                     20  q.rx  Minimum range factor per iteration [0.5]
          xp{n}(:) list of available probe SNR values
          lf       log file ID
          i        model probed
          x        probe SNR value used
          r        response: 0=wrong, 1=right.
          o        plot option string:
                   'p' Plot pdf
                   'h' Plot probe history
                   'x' Plot expected cost function for probe SNRs
                   'c' Plot cost function evolution

 Outputs:
          xx       recommended probe SNR
          ii       recommended model to probe next
          m(2,n,3) estimated srt and slope of all models
                   m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
          v(3,n)   estimated covariance matrix entries:
                   [var(srt) cov(srt,slope) var(slope)]'
          mr(2,n,3)robust estimated srt and slope of all models
                   m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
          vr(3,n)  robust estimated covariance matrix entries:
                   [var(srt) cov(srt,slope) var(slope)]'
          msr(:,3)  List probe snrs and results: [model probe-snr result]

 Algorithm parameters:

 The algorithm estimates the psychometric function for a number of models simultaneously.
 A psychometric function specifies the probability of correct recognition as a function of
 SNR and is completely specified by two free parameters: the SNR (in dB) and the slope (in 1/dB)
 at a specified target recognition probability (e.g. 0.5 or 0.75). The p(:,n) parameters specify
 some of the details of the psychometric function and can be different for each model:
   p(1,n) gives the target recognition probability
   p(2,n) gives the probability of incorrect recognition at very good SNRs (the &quot;miss&quot; probability).
          If this value is made too low, a single unwarrented recognition error will have a large
          effect of the estimated parameters and greatly lengthen the time to convergence. The default
          value is 4%.
   p(3,n) gives the probability of correct recognition at very poor SNRs (the &quot;guess&quot; probabiity).
          This should be set to 1/N where N is the number of possible responses to a stimulus.
 p(4:5,n) These give the initial min and max SNR at the target recognition probability. They will
           be adjusted by the algorithm if necessary but wrong values will delay convergence.
 p(6:7,n) These given the initial min and max slope (in probability per dB) at the target
          recognition probability.
 The remaining parameters are shared between all the models and control how the algorithm operates.
   q(1:2) These parameters determine the sampling size of the joint pdf in SNR and Slope respectively.
   q(3)   This parameter specifies how many potential values of probe SNR to evaluate in order to
          determine which is likely to give the most improvement in the parameter estimates.
   q(4)   This parameter gives the relative weight of SNR and Slope in the cost function. Increasing
          its value will improve the slope estimate (or log(slope) estimate) at the expense of the
          SNR estimate. Actually its value is not all that critical.
   q(5)   At each iteration v_psycest evaluates several potential probe SNR levels to see which is
          likely to give the most improvement to the parameter estimates. This parameter fixes the
          minimum spacing between these potential probe values. It only has an effect when the variances
          of the parameter estimates are very small.
   q(6)   To prevent the routine getting lost, this parameter specifies the smallest reasonable value
          of the Slope parameter at the target recognition probability. Under normal circumstances, it
          has no effect.
   q(7)   For each model, the routine maintains a sampled version of the joint pdf of the SNR and slope.
          The sampling grid is reclculated after each iteration to encompass this number of standard
          deviations in each dimension.
   q(8)   At each iteration, v_psycest advises which model to probe next on the basis of which
          gives the greatest expected reduction in cost function. To ensure that all models
          are periodically sampled, this expected reduction of each unprobed model is multiplied
          by this parameter at each iteration. Thus a model will eventually be probed even if
          its expected cost factor improvement is small.
   q(9)   This determines whether the recommended probe SNR to use at the next iteration is chosen to
          minimize the expected variance or the v_entropy of the estimated SNR and slope distributions.
          My experience is that v_entropy (the default) gives faster convergence.
  q(10)   This selects whether the underlying model is a logistic function (1) or a cumulative
          gaussian (2). The choice makes little difference unless the &quot;miss&quot; or &quot;guess&quot; probabilities
          are very very small.
  q(11)   This selects whether the internal pdf samples &quot;slope&quot; (0) or &quot;log(slope)&quot; (1). It is recommended
          that you stick to the default of log(slope) since this results in more symmetrical distributions. [1]
  q(12)   This sets the number of std deviations of the prior corresponding to the limits of the
          initial distribution set in p(4:7,n). A small value corresponds to a very weak prior. [1]
  q(13)   Probability floor (integrated over entire grid) [0.0001]
  q(14)   Number of std devs to explore when evaluating possible probe SNRs [2]
  q(15)   Maximum probe SNR shift beyond the range of previous probes (dB) [10]
  q(16)   Grid interpolation threshold. The grid is only interpolated if
          it has shrunk by less that this factor [0.5]
  q(17)   If the desired grid limits change by less than this fraction of
          the grid length, they will be left unaltered. [0.1]
  q(18)   Look 2-ahead when choosing probe [1]
  q(19)   Outlier probability; probe results with lower probability than
          this will be ignored [0.01]
  q(20)   The ranges of SRT and log slope will not be reduced per iteration below this factor [0.5]</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../matlabicon.gif)">
<li><a href="v_axisenlarge.html" class="code" title="function v_axisenlarge(f,h)">v_axisenlarge</a>	V_AXISENLARGE - enlarge the axes of a figure (f,h)</li><li><a href="v_quadpeak.html" class="code" title="function [v,x,t,m,ze]=v_quadpeak(z)">v_quadpeak</a>	V_PEAK2DQUAD find quadratically-interpolated peak in a N-D array</li><li><a href="v_texthvc.html" class="code" title="function h=v_texthvc(x,y,t,p,q,r)">v_texthvc</a>	V_TEXTHVC - write text on graph with specified alignment and colour</li></ul>
This function is called by:
<ul style="list-style-image:url(../matlabicon.gif)">
<li><a href="v_psycdigit.html" class="code" title="function [m,v]=v_psycdigit(proc,r,mode,p,q,xp,noise,fn,dfile,ofile)">v_psycdigit</a>	V_PSYCDIGIT measures psychometric function using TIDIGITS stimuli</li></ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [xx,ii,m,v,mr,vr]=v_psycest(iq,x,r,xp,lf)</a>
0002 <span class="comment">%V_PSYCEST estimate multiple psychometric functions</span>
0003 <span class="comment">%</span>
0004 <span class="comment">% Usage: [xx,ii,m,v]=v_psycest(-n,p,q,xp,lf) % initialize n models</span>
0005 <span class="comment">%        [xx,ii,m,v]=v_psycest(i,x,r)     % supply a trial result to v_psycest</span>
0006 <span class="comment">%        [xx,ii,m,v]=v_psycest(i,x,r,o)   % supply a trial result to v_psycest and plot</span>
0007 <span class="comment">%                    v_psycest(i,o)       % plot pdf of model i with plot options o</span>
0008 <span class="comment">%        [xx,ii,m,v,mr,vr]=v_psycest(i,o) % Determine robust outputs in addition expcluding outliers [takes longer]</span>
0009 <span class="comment">%              [p,q,msr]=psychest(0)    % output model parameters (or print them if no outputs)</span>
0010 <span class="comment">%</span>
0011 <span class="comment">% Inputs (see usage examples for argument positions):</span>
0012 <span class="comment">%         -n        minus the number of models</span>
0013 <span class="comment">%          p(:,n)   parameters for each model</span>
0014 <span class="comment">%                      1  thresh [0.75]</span>
0015 <span class="comment">%                      2  miss [0.04]</span>
0016 <span class="comment">%                      3  guess [0.1]</span>
0017 <span class="comment">%                      4  SNR min [-20]</span>
0018 <span class="comment">%                      5  SNR max [20]</span>
0019 <span class="comment">%                      6  Slope min (but over-ridden if &lt;q.sl) [0]</span>
0020 <span class="comment">%                      7  slope max [0.5]</span>
0021 <span class="comment">%          q(:)     parameters common to all models (vector or struct)</span>
0022 <span class="comment">%                      1  q.nx  number of SNR values in pdf [40]</span>
0023 <span class="comment">%                      2  q.ns  number of slope values in pdf [21]</span>
0024 <span class="comment">%                      3  q.nh  number of probe SNR values to evaluate [30]</span>
0025 <span class="comment">%                      4  q.cs  weighting of slope relative to SRT in cost function [1]</span>
0026 <span class="comment">%                      5  q.dh  minimum step size in dB for probe SNRs [0.2]</span>
0027 <span class="comment">%                      6  q.sl  min slope at threshold (must be &gt;0) [0.005]</span>
0028 <span class="comment">%                      7  q.kp  number of std deviations of the pdfs to keep [4]</span>
0029 <span class="comment">%                      8  q.hg  amount to grow expected gains in ni trials [1.3]</span>
0030 <span class="comment">%                      9  q.cf  cost function: 1=variance, 2=v_entropy [2]</span>
0031 <span class="comment">%                     10  q.pm  psychometric model: 1=logistic, 2=cumulative gaussian [1]</span>
0032 <span class="comment">%                     11  q.lg  use log slope in pdf: 0=no, 1=yes [1]</span>
0033 <span class="comment">%                     12  q.pp  Number of prior standard deviations in initial semi-range [1]</span>
0034 <span class="comment">%                     13  q.pf  Probability floor (integrated over entire grid) [0.0001]</span>
0035 <span class="comment">%                     14  q.ts  Number of std devs to explore [2]</span>
0036 <span class="comment">%                     15  q.dp  Maximum probe SNR shift (dB) [10]</span>
0037 <span class="comment">%                     16  q.it  Grid interpolation threshold [0.5]</span>
0038 <span class="comment">%                     17  q.at  Axis change threshold [0.1]</span>
0039 <span class="comment">%                     18  q.la  Look 2-ahead when choosing probe [1]</span>
0040 <span class="comment">%                     19  q.op  Outlier probability [0.01]</span>
0041 <span class="comment">%                     20  q.rx  Minimum range factor per iteration [0.5]</span>
0042 <span class="comment">%          xp{n}(:) list of available probe SNR values</span>
0043 <span class="comment">%          lf       log file ID</span>
0044 <span class="comment">%          i        model probed</span>
0045 <span class="comment">%          x        probe SNR value used</span>
0046 <span class="comment">%          r        response: 0=wrong, 1=right.</span>
0047 <span class="comment">%          o        plot option string:</span>
0048 <span class="comment">%                   'p' Plot pdf</span>
0049 <span class="comment">%                   'h' Plot probe history</span>
0050 <span class="comment">%                   'x' Plot expected cost function for probe SNRs</span>
0051 <span class="comment">%                   'c' Plot cost function evolution</span>
0052 <span class="comment">%</span>
0053 <span class="comment">% Outputs:</span>
0054 <span class="comment">%          xx       recommended probe SNR</span>
0055 <span class="comment">%          ii       recommended model to probe next</span>
0056 <span class="comment">%          m(2,n,3) estimated srt and slope of all models</span>
0057 <span class="comment">%                   m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates</span>
0058 <span class="comment">%          v(3,n)   estimated covariance matrix entries:</span>
0059 <span class="comment">%                   [var(srt) cov(srt,slope) var(slope)]'</span>
0060 <span class="comment">%          mr(2,n,3)robust estimated srt and slope of all models</span>
0061 <span class="comment">%                   m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates</span>
0062 <span class="comment">%          vr(3,n)  robust estimated covariance matrix entries:</span>
0063 <span class="comment">%                   [var(srt) cov(srt,slope) var(slope)]'</span>
0064 <span class="comment">%          msr(:,3)  List probe snrs and results: [model probe-snr result]</span>
0065 <span class="comment">%</span>
0066 <span class="comment">% Algorithm parameters:</span>
0067 <span class="comment">%</span>
0068 <span class="comment">% The algorithm estimates the psychometric function for a number of models simultaneously.</span>
0069 <span class="comment">% A psychometric function specifies the probability of correct recognition as a function of</span>
0070 <span class="comment">% SNR and is completely specified by two free parameters: the SNR (in dB) and the slope (in 1/dB)</span>
0071 <span class="comment">% at a specified target recognition probability (e.g. 0.5 or 0.75). The p(:,n) parameters specify</span>
0072 <span class="comment">% some of the details of the psychometric function and can be different for each model:</span>
0073 <span class="comment">%   p(1,n) gives the target recognition probability</span>
0074 <span class="comment">%   p(2,n) gives the probability of incorrect recognition at very good SNRs (the &quot;miss&quot; probability).</span>
0075 <span class="comment">%          If this value is made too low, a single unwarrented recognition error will have a large</span>
0076 <span class="comment">%          effect of the estimated parameters and greatly lengthen the time to convergence. The default</span>
0077 <span class="comment">%          value is 4%.</span>
0078 <span class="comment">%   p(3,n) gives the probability of correct recognition at very poor SNRs (the &quot;guess&quot; probabiity).</span>
0079 <span class="comment">%          This should be set to 1/N where N is the number of possible responses to a stimulus.</span>
0080 <span class="comment">% p(4:5,n) These give the initial min and max SNR at the target recognition probability. They will</span>
0081 <span class="comment">%           be adjusted by the algorithm if necessary but wrong values will delay convergence.</span>
0082 <span class="comment">% p(6:7,n) These given the initial min and max slope (in probability per dB) at the target</span>
0083 <span class="comment">%          recognition probability.</span>
0084 <span class="comment">% The remaining parameters are shared between all the models and control how the algorithm operates.</span>
0085 <span class="comment">%   q(1:2) These parameters determine the sampling size of the joint pdf in SNR and Slope respectively.</span>
0086 <span class="comment">%   q(3)   This parameter specifies how many potential values of probe SNR to evaluate in order to</span>
0087 <span class="comment">%          determine which is likely to give the most improvement in the parameter estimates.</span>
0088 <span class="comment">%   q(4)   This parameter gives the relative weight of SNR and Slope in the cost function. Increasing</span>
0089 <span class="comment">%          its value will improve the slope estimate (or log(slope) estimate) at the expense of the</span>
0090 <span class="comment">%          SNR estimate. Actually its value is not all that critical.</span>
0091 <span class="comment">%   q(5)   At each iteration v_psycest evaluates several potential probe SNR levels to see which is</span>
0092 <span class="comment">%          likely to give the most improvement to the parameter estimates. This parameter fixes the</span>
0093 <span class="comment">%          minimum spacing between these potential probe values. It only has an effect when the variances</span>
0094 <span class="comment">%          of the parameter estimates are very small.</span>
0095 <span class="comment">%   q(6)   To prevent the routine getting lost, this parameter specifies the smallest reasonable value</span>
0096 <span class="comment">%          of the Slope parameter at the target recognition probability. Under normal circumstances, it</span>
0097 <span class="comment">%          has no effect.</span>
0098 <span class="comment">%   q(7)   For each model, the routine maintains a sampled version of the joint pdf of the SNR and slope.</span>
0099 <span class="comment">%          The sampling grid is reclculated after each iteration to encompass this number of standard</span>
0100 <span class="comment">%          deviations in each dimension.</span>
0101 <span class="comment">%   q(8)   At each iteration, v_psycest advises which model to probe next on the basis of which</span>
0102 <span class="comment">%          gives the greatest expected reduction in cost function. To ensure that all models</span>
0103 <span class="comment">%          are periodically sampled, this expected reduction of each unprobed model is multiplied</span>
0104 <span class="comment">%          by this parameter at each iteration. Thus a model will eventually be probed even if</span>
0105 <span class="comment">%          its expected cost factor improvement is small.</span>
0106 <span class="comment">%   q(9)   This determines whether the recommended probe SNR to use at the next iteration is chosen to</span>
0107 <span class="comment">%          minimize the expected variance or the v_entropy of the estimated SNR and slope distributions.</span>
0108 <span class="comment">%          My experience is that v_entropy (the default) gives faster convergence.</span>
0109 <span class="comment">%  q(10)   This selects whether the underlying model is a logistic function (1) or a cumulative</span>
0110 <span class="comment">%          gaussian (2). The choice makes little difference unless the &quot;miss&quot; or &quot;guess&quot; probabilities</span>
0111 <span class="comment">%          are very very small.</span>
0112 <span class="comment">%  q(11)   This selects whether the internal pdf samples &quot;slope&quot; (0) or &quot;log(slope)&quot; (1). It is recommended</span>
0113 <span class="comment">%          that you stick to the default of log(slope) since this results in more symmetrical distributions. [1]</span>
0114 <span class="comment">%  q(12)   This sets the number of std deviations of the prior corresponding to the limits of the</span>
0115 <span class="comment">%          initial distribution set in p(4:7,n). A small value corresponds to a very weak prior. [1]</span>
0116 <span class="comment">%  q(13)   Probability floor (integrated over entire grid) [0.0001]</span>
0117 <span class="comment">%  q(14)   Number of std devs to explore when evaluating possible probe SNRs [2]</span>
0118 <span class="comment">%  q(15)   Maximum probe SNR shift beyond the range of previous probes (dB) [10]</span>
0119 <span class="comment">%  q(16)   Grid interpolation threshold. The grid is only interpolated if</span>
0120 <span class="comment">%          it has shrunk by less that this factor [0.5]</span>
0121 <span class="comment">%  q(17)   If the desired grid limits change by less than this fraction of</span>
0122 <span class="comment">%          the grid length, they will be left unaltered. [0.1]</span>
0123 <span class="comment">%  q(18)   Look 2-ahead when choosing probe [1]</span>
0124 <span class="comment">%  q(19)   Outlier probability; probe results with lower probability than</span>
0125 <span class="comment">%          this will be ignored [0.01]</span>
0126 <span class="comment">%  q(20)   The ranges of SRT and log slope will not be reduced per iteration below this factor [0.5]</span>
0127 
0128 <span class="comment">%      Copyright (C) Mike Brookes and Clement Doire 2009-2016</span>
0129 <span class="comment">%      Version: $Id: v_psycest.m 10865 2018-09-21 17:22:45Z dmb $</span>
0130 <span class="comment">%</span>
0131 <span class="comment">%   VOICEBOX is a MATLAB toolbox for speech processing.</span>
0132 <span class="comment">%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html</span>
0133 <span class="comment">%</span>
0134 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0135 <span class="comment">%   This program is free software; you can redistribute it and/or modify</span>
0136 <span class="comment">%   it under the terms of the GNU General Public License as published by</span>
0137 <span class="comment">%   the Free Software Foundation; either version 2 of the License, or</span>
0138 <span class="comment">%   (at your option) any later version.</span>
0139 <span class="comment">%</span>
0140 <span class="comment">%   This program is distributed in the hope that it will be useful,</span>
0141 <span class="comment">%   but WITHOUT ANY WARRANTY; without even the implied warranty of</span>
0142 <span class="comment">%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the</span>
0143 <span class="comment">%   GNU General Public License for more details.</span>
0144 <span class="comment">%</span>
0145 <span class="comment">%   You can obtain a copy of the GNU General Public License from</span>
0146 <span class="comment">%   http://www.gnu.org/copyleft/gpl.html or by writing to</span>
0147 <span class="comment">%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.</span>
0148 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0149 
0150 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0151 <span class="comment">% Bugs/Suggestions:</span>
0152 <span class="comment">% (2)use a structure for input parameters</span>
0153 <span class="comment">% (3) v_entropy seems to change signifiacntly when the grid is rescaled</span>
0154 <span class="comment">% (4)add a forgetting factor for old measurements</span>
0155 <span class="comment">% (8) use quadratic interpolation to find cost function minimum</span>
0156 <span class="comment">% (10) optionally output the whole pdf + its axis values</span>
0157 <span class="comment">% (11) optionally output all model probe values</span>
0158 <span class="comment">% (13) Should perhaps estimate and return the mean and slope compensated for the guess rate</span>
0159 <span class="comment">%      i.e. if you want p, you actually test at guess+ p*(1-guess-miss)/(1-miss)</span>
0160 <span class="comment">% (17) remember probe snrs, models and results and recalculate the entire</span>
0161 <span class="comment">%      array when changing precision</span>
0162 <span class="comment">% (18) could force probe SNRs to be a multiple of minimum step size</span>
0163 <span class="comment">% (19) use a non-uniform prior e.e. 1/(1+x^2)</span>
0164 <span class="comment">% (20) possibly use different parametrization (e.g. 1/slope or a second SRT threshold)</span>
0165 <span class="comment">% (22) save inputs so that pdfs can be recalculated when rescaling</span>
0166 <span class="comment">% (24) Parametrize slope as x=(100s^2-1)/20s to make the distribution more uniform</span>
0167 <span class="comment">%      inverse is s=(sqrt(x^2+1)-x)/10; alternatively use log(slope)</span>
0168 <span class="comment">% (25) optionally have no prior (i.e. maximum likelihood)</span>
0169 <span class="comment">% (26) Check why the estimated variances are too small</span>
0170 <span class="comment">% (27) Adapt range of potential probes if optimum is at the limit</span>
0171 <span class="comment">% (28) Resample the pdf on the basis of the marginal distributions</span>
0172 <span class="comment">% (29) Select the probe range on the basis of marginal distributions</span>
0173 <span class="comment">% (30) Resample the pdfs by a factor of 2 always</span>
0174 <span class="comment">% (31) use uneven samples in pdf concentrated in the middle</span>
0175 <span class="comment">% (32) Use a parametric mixture including one-sided constants for the pdf</span>
0176 <span class="comment">% (33) Selection of subset of prescribed probe SNRs is not optimum</span>
0177 <span class="comment">% (34) use quadratic interpolation to select the probe SNR unless using fixed values</span>
0178 <span class="comment">% (35) expand the pdf grid based on the effective number of samples (like particle filters)</span>
0179 
0180 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0181 <span class="comment">% persistent variables</span>
0182 <span class="comment">% sizes: nx=SNR values in pdf</span>
0183 <span class="comment">%        ns=slope values in pdf (1/std dev of cumulative gaussian)</span>
0184 <span class="comment">%        nh=potential test SNRs</span>
0185 <span class="comment">%        ni=simultaneous models being estimated</span>
0186 <span class="comment">% wq(ns*nx,ni)  = log prob of model; each column is a vectorized ns*nx matrix</span>
0187 <span class="comment">% nr(1,10)      = parameter values: [SNR-pdf slope-pdf SNR-probe ...]</span>
0188 <span class="comment">% pr(7,ni)      = input model parameters</span>
0189 <span class="comment">% qr(4,ni)      = derived model parameters</span>
0190 <span class="comment">% xq(nr(1),ni)  = SNR values in pdf</span>
0191 <span class="comment">% sq(nr(2),ni)  = slope values in pdf</span>
0192 <span class="comment">% mq(2,ni,3)    = estimated srt and slope of all models</span>
0193 <span class="comment">% vq(3,ni)      = estimated covariance matrix entries:[var(srt) cov(srt,slope) var(slope)]'</span>
0194 <span class="comment">% xn(1,ni)      = next probe value to use</span>
0195 <span class="comment">% hn(1,ni)      = expected decrease in cost function after next probe</span>
0196 <span class="comment">% hfact</span>
0197 <span class="comment">% xz</span>
0198 <span class="comment">% res(nres,7)   = results [model-number, probe-SNR, result]</span>
0199 <span class="comment">% nres          = total number of probe results</span>
0200 <span class="comment">% nresq(1,ni)   = number of probe results for each model</span>
0201 <span class="comment">% xmm(2,ni)     = [min(probe-SNR); max(probe-SNR)]</span>
0202 <span class="comment">% mq0(2,ni)      = prior means [x-mean; s-mean]</span>
0203 <span class="comment">% pq0(2,ni)      = scale factors for prior distribution</span>
0204 <span class="comment">% wfl           = floor relative to peak of log-pdf array</span>
0205 <span class="comment">% sqmin         = minimum value of slope or log-slope</span>
0206 <span class="comment">% LOG           = file ID of logging file</span>
0207 <span class="comment">% mqr(2,ni,3)   = robust mean estimates (excluding outliers)</span>
0208 <span class="comment">% vqr(3,ni)     = robust variance estimates (excluding outliers)</span>
0209 <span class="comment">% nresr         = nres value at which robust estimates were last calculated</span>
0210 
0211 <span class="keyword">persistent</span> wq xq sq nr pr qr mq vq xn hn hfact xz res nres nresq xmm mq0 pq0 wfl sqmin LOG mqr vqr nresr xlim
0212 
0213 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0214 <span class="comment">%              Initialization (iq&lt;0)            %</span>
0215 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0216 <span class="keyword">if</span> iq&lt;0  <span class="comment">% initialization</span>
0217     <span class="keyword">if</span> nargin&gt;4 &amp;&amp; ~isempty(lf)
0218         LOG=lf;
0219         fprintf(LOG,<span class="string">'******************************\npsycest Initialization: %s\n******************************\n'</span>,datestr(now));
0220     <span class="keyword">else</span>
0221         LOG=0;
0222     <span class="keyword">end</span>
0223     ni=-iq;                 <span class="comment">% number of models</span>
0224     nres=0;                 <span class="comment">% total number of probe-results so far</span>
0225     nresq=zeros(1,ni);      <span class="comment">% number of probe-results for each model</span>
0226     res=zeros(1,7);         <span class="comment">% array for saving results [expands as needed]</span>
0227     pr=repmat([0.75 0.04 0.1 -20 20 0 0.5]',1,ni);          <span class="comment">% default parameters</span>
0228     <span class="keyword">if</span> nargin&gt;1
0229         <span class="keyword">if</span> size(x,2)&gt;1
0230             <span class="keyword">if</span> size(x,2)&gt;ni
0231                 error(<span class="string">'initialization parameter argument has too many columns'</span>);
0232             <span class="keyword">end</span>
0233             pr(1:size(x,1),1:size(x,2))=x;
0234         <span class="keyword">else</span>
0235             pr(1:size(x,1),:)=repmat(x,1,ni);
0236         <span class="keyword">end</span>
0237     <span class="keyword">end</span>
0238     nr=[40 21 30 1 0.2 0.02 4 1.3 2 1 1 1 0.0001 2 10 0.5 0.1 1 0.01 0.5]';                      <span class="comment">% default parameter values</span>
0239     nrf={<span class="string">'nx'</span>,<span class="string">'ns'</span>,<span class="string">'nh'</span>,<span class="string">'cs'</span>,<span class="string">'dh'</span>,<span class="string">'sl'</span>,<span class="string">'kp'</span>,<span class="string">'hg'</span>,<span class="string">'cf'</span>,<span class="string">'pm'</span>,<span class="string">'lg'</span>,<span class="string">'pp'</span>, <span class="string">'pf'</span>, <span class="string">'ts'</span>, <span class="string">'dp'</span>, <span class="string">'it'</span>,<span class="string">'at'</span>,<span class="string">'la'</span>,<span class="string">'op'</span>,<span class="string">'rx'</span>};     <span class="comment">% parameter field names</span>
0240     numnr=length(nr);
0241     <span class="keyword">if</span> nargin&gt;2 &amp;&amp; ~isempty(r) <span class="comment">% replace defaults with any parameter values specified in the call</span>
0242         <span class="keyword">if</span> isstruct(r)
0243             fnn=fieldnames(r);
0244             <span class="keyword">for</span> i=1:length(fnn)
0245                 mk=strcmp(fnn{i},nrf);
0246                 <span class="keyword">if</span> any(mk)
0247                     nr(mk)=r.(fnn{i});
0248                 <span class="keyword">end</span>
0249             <span class="keyword">end</span>
0250         <span class="keyword">else</span>
0251             nr(1:min(numnr,numel(r)))=r(:); <span class="comment">% if parameters are specified as a vector, copy them across</span>
0252         <span class="keyword">end</span>
0253         nr(1:3)=round(nr(1:3));             <span class="comment">% first three parameters must be integers</span>
0254     <span class="keyword">end</span>
0255     pr(6,:)=max(pr(6,:),nr(6));                <span class="comment">% low limit of slope in prob/dB</span>
0256     nxq=nr(1);
0257     nsq=nr(2);
0258     nsxq=nxq*nsq;
0259     xq=(0:nxq-1)'*(pr(5,:)-pr(4,:))/(nxq-1)+repmat(pr(4,:),nxq,1);
0260     <span class="keyword">if</span> nr(11)  <span class="comment">% if log slope</span>
0261         sqmin=log(nr(6)); <span class="comment">% low limit for axis</span>
0262         sq=(1-nsq:0)'*(log(pr(7,:))-log(max(pr(6,:),nr(6))))/(nsq-1)+repmat(log(pr(7,:)),nsq,1);
0263     <span class="keyword">else</span>
0264         sqmin=nr(6); <span class="comment">% low limit for axis</span>
0265         sq=(1-nsq:0)'*(pr(7,:)-max(pr(6,:),nr(6)))/(nsq-1)+repmat(pr(7,:),nsq,1);
0266     <span class="keyword">end</span>
0267     mq0=[mean(xq,1); mean(sq,1)];     <span class="comment">% find mean of prior</span>
0268     pq0=-2*nr(12)^2*[xq(<span class="keyword">end</span>,:)-xq(1,:); sq(<span class="keyword">end</span>,:)-sq(1,:)].^(-2); <span class="comment">% scale factor for prior distribution</span>
0269     wq=repmat(pq0(2,1)*(sq(:,1)-mq0(2,1)).^2,1,nxq)+repmat(pq0(1,1)*(xq(:,1)-mq0(1,1))'.^2,nsq,1); <span class="comment">% log probs of prior (same for all models)</span>
0270     wfl=log(nr(13)/nsxq); <span class="comment">% floor relative to peak of wq array</span>
0271     wq=repmat(wq(:),1,ni);              <span class="comment">% initialize to +-q.pp std deviations of gaussian prior</span>
0272     qr=zeros(5,ni);
0273     qr(1,:)=1-pr(2,:)-pr(3,:);             <span class="comment">% prob range covered by cumulative gaussian</span>
0274     <span class="keyword">if</span> any(qr(1,:)&lt;=0)
0275         [dum,i]=min(qr(1,:));
0276         error(<span class="string">'Model %d: guess prob (%.2f) + miss prob (%.2f) is &gt;=1'</span>,i,pr(3,i),pr(2,:));
0277     <span class="keyword">end</span>
0278     qr(2,:)=(pr(1,:)-pr(3,:))./qr(1,:);                 <span class="comment">% cumulative gaussian probability at threshold</span>
0279     <span class="keyword">if</span> any(abs(qr(2,:)-0.5)&gt;=0.5)
0280         [dum,i]=max(qr(2,:)-0.5);
0281         error(<span class="string">'Model %d: target SRT threshold (%.2f) must lie in the range guess prob (%.2f) to (1-miss prob) (%.2f)'</span>,i,pr(1,i),pr(3,i),1-pr(2,:));
0282     <span class="keyword">end</span>
0283     <span class="keyword">switch</span> nr(10)                                       <span class="comment">% switch on model type</span>
0284         <span class="keyword">case</span> 1                                          <span class="comment">% logistic model</span>
0285             qr(3,:)=log(qr(2,:)./(1-qr(2,:)));            <span class="comment">% x position of target in std measure</span>
0286             qr(4,:)=qr(1,:).*qr(2,:).*(1-qr(2,:));        <span class="comment">% slope*&quot;stddev&quot; at threshold</span>
0287         <span class="keyword">case</span> 2                                          <span class="comment">% cumulative gaussian model</span>
0288             qr(3,:)=norminv(qr(2,:));                   <span class="comment">% x position of target in std measure</span>
0289             qr(4,:)=qr(1,:).*normpdf(qr(3,:));           <span class="comment">% slope*stddev at threshold</span>
0290         <span class="keyword">otherwise</span>
0291             error(<span class="string">'Unrecognised psychometric model selection'</span>);
0292     <span class="keyword">end</span>
0293     <span class="keyword">if</span> nr(9)&lt; 1 || nr(9)&gt;2
0294         error(<span class="string">'Unrecognised cost function option'</span>);
0295     <span class="keyword">end</span>
0296     mq=repmat(mq0,[1,1,3]);        <span class="comment">% initial means, joint mode and marginal mode all are equal</span>
0297     vq=[var(xq,1,1); zeros(1,ni); var(sq,1,1)];         <span class="comment">% initial variances (actually ignores prior probs)</span>
0298     mqr=mq;         <span class="comment">% robust means and modes</span>
0299     vqr=vq;         <span class="comment">% robust variances</span>
0300     nresr=0;        <span class="comment">% robust calculation time</span>
0301     xlim=[-Inf(1,ni); Inf(1,ni)]; <span class="comment">% SNR limits for outliers</span>
0302     <span class="comment">%     hn=[1-nr(4) 0 nr(4)]*vq;  % artificially high value of cost function ensures all models are probed early on</span>
0303     hn=Inf(1,ni);           <span class="comment">% very high initial cost function</span>
0304     hfact=nr(8)^(1/ni);     <span class="comment">% growth factor to ensure that no model is neglected for too long</span>
0305     xn=mq0(1,:);             <span class="comment">% start at mean value</span>
0306     xmm=[xn; xn];  <span class="comment">% min/max of probe values for each model</span>
0307     <span class="keyword">if</span> nargin&gt;=4 &amp;&amp; ~isempty(xp)
0308         <span class="keyword">if</span> iscell(xp)
0309             xz=xp;
0310         <span class="keyword">else</span>
0311             xz=repmat(num2cell(xp(:)',2),ni,1);        <span class="comment">% else replicate for each model</span>
0312         <span class="keyword">end</span>
0313         <span class="keyword">for</span> i=1:ni
0314             [dum,j]=min(abs(xz{i}-mq0(1,i)));      <span class="comment">% select the closest available probe to the mean</span>
0315             xn(i)=xz{i}(j);
0316         <span class="keyword">end</span>
0317     <span class="keyword">else</span>
0318         xz=cell(ni,1);          <span class="comment">% make empty cells</span>
0319     <span class="keyword">end</span>
0320 <span class="keyword">elseif</span> iq&gt;0
0321     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0322     <span class="comment">%     iq&gt;0 means update or plot model iq           %</span>
0323     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0324     <span class="comment">%        first check the plotting options          %</span>
0325     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0326     <span class="keyword">if</span> nargin&lt;3
0327         <span class="keyword">if</span> nargin==2
0328             po=x; <span class="comment">% plot options are 2nd argument</span>
0329         <span class="keyword">elseif</span> ~nargout
0330             po=<span class="string">'p'</span>;
0331         <span class="keyword">else</span>
0332             po=<span class="string">''</span>;
0333         <span class="keyword">end</span>
0334     <span class="keyword">else</span>
0335         <span class="keyword">if</span> nargin&gt;3
0336             po=xp; <span class="comment">% plot options are 4th argument</span>
0337         <span class="keyword">else</span>
0338             po=<span class="string">''</span>;
0339         <span class="keyword">end</span>
0340     <span class="keyword">end</span>
0341     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0342     <span class="comment">%     now get parameters of current model (iq)     %</span>
0343     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0344     nxq=nr(1);                  <span class="comment">% number of SNR values</span>
0345     nsq=nr(2);                  <span class="comment">% number of slope values</span>
0346     nxh=nr(3);                  <span class="comment">% number of probe-SNR values</span>
0347     nsxq=nxq*nsq;               <span class="comment">% size of log pdf array</span>
0348     thresh=pr(1,iq);            <span class="comment">% target probability at threshold</span>
0349     guess=pr(3,iq);             <span class="comment">% guess rate (1/choices)</span>
0350     pscale=qr(1,iq);            <span class="comment">% prob range left after subtracting miss and guess probs</span>
0351     xtstd=qr(3,iq);             <span class="comment">% x position of target in std measure</span>
0352     xqi=xq(:,iq);               <span class="comment">% SRT values of the pdf array</span>
0353     sqi=sq(:,iq);               <span class="comment">% slope (or log slope) values in PDF</span>
0354     wqi=wq(:,iq);                  <span class="comment">% log pdf array</span>
0355     mqi= mq(:,iq,1);            <span class="comment">% [xe; se] means</span>
0356     vqi=vq(:,iq);               <span class="comment">% [xv; sxv; sv]  covariance matrix</span>
0357     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0358     <span class="comment">%     rescale the pdfs if necessary         %</span>
0359     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0360     ssd=sqrt(vqi(3));           <span class="comment">% std deviation of slope</span>
0361     xsd=sqrt(vqi(1));           <span class="comment">% std deviation of SRT</span>
0362     <span class="comment">% determine range of new grid</span>
0363     xqrange=xqi(end)-xqi(1);    <span class="comment">% range of old SRT grid</span>
0364     sqrange=sqi(end)-sqi(1);    <span class="comment">% range of old slope grid</span>
0365     <span class="keyword">if</span>  nresq(iq)&lt;2              <span class="comment">% keep old limits if nresq(iq)&lt;2</span>
0366         xq2lim=xqi([1 end]);
0367         sq2lim=sqi([1 end]);
0368     <span class="keyword">else</span>
0369         xqsemirange=max(nr(7)*xsd,0.5*nr(20)*xqrange);
0370         xq2lim=[mqi(1)-xqsemirange, mqi(1)+xqsemirange];
0371         sqsemirange=max(nr(7)*ssd,0.5*nr(20)*sqrange);
0372         sq2lim=[max(sqmin,mqi(2)-sqsemirange),mqi(2)+sqsemirange];
0373         <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0374         <span class="comment">% To avoid unnecessary changes, keep old limits if the change  %</span>
0375         <span class="comment">% is less than nr(17) times the previous range [nr(17)=0.1]    %</span>
0376         <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0377         <span class="keyword">if</span> abs(xq2lim(1)-xqi(1))&lt;nr(17)*xqrange
0378             xq2lim(1)=xqi(1);
0379         <span class="keyword">end</span>
0380         <span class="keyword">if</span> abs(xq2lim(2)-xqi(end))&lt;nr(17)*xqrange
0381             xq2lim(2)=xqi(end);
0382         <span class="keyword">end</span>
0383         
0384         <span class="keyword">if</span> abs(sq2lim(1)-sqi(1))&lt;nr(17)*sqrange
0385             sq2lim(1)=sqi(1);
0386         <span class="keyword">end</span>
0387         <span class="keyword">if</span> abs(sq2lim(2)-sqi(end))&lt;nr(17)*sqrange
0388             sq2lim(2)=sqi(end);
0389         <span class="keyword">end</span>
0390     <span class="keyword">end</span>
0391     xq2=linspace(xq2lim(1),xq2lim(2),nxq)';   <span class="comment">% new x axis values</span>
0392     sq2=linspace(sq2lim(1),sq2lim(2),nsq)';   <span class="comment">% new s axis values</span>
0393     wqup=2;                                   <span class="comment">% update flag</span>
0394     <span class="keyword">if</span> xq2(1)&lt;xqi(1) || xq2(end)&gt;xqi(end) || sq2(1)&lt;sqi(1) || sq2(end)&gt;sqi(end)
0395         <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0396         <span class="comment">% if extrapolating, recalculate log-pdfs from saved data  %</span>
0397         <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0398         <span class="keyword">if</span> LOG, fprintf(LOG,<span class="string">'N=%d:%d, extrapolate: x(%.2g,%.2g)-&gt;(%.2g,%.2g) and s(%.2g,%.2g)-&gt;(%.2g,%.2g)\n'</span>,iq,nresq(iq),xqi([1 end]),xq2([1 end]),sqi([1 end]),sq2([1 end])); <span class="keyword">end</span>
0399         wq2=repmat(pq0(2,iq)*(sq2-mq0(2,iq)).^2,1,nxq)+repmat(pq0(1,iq)*(xq2-mq0(1,iq))'.^2,nsq,1); <span class="comment">% log probs of prior for model iq</span>
0400         ires=find(res(1:nres,1)==iq); <span class="comment">% find the results for this model</span>
0401         <span class="keyword">if</span> nr(11)                   <span class="comment">% if log slope</span>
0402             sqis=exp(sq2)/qr(4,iq);    <span class="comment">% inverse std dev of gaussian (proportional to slope)</span>
0403         <span class="keyword">else</span>
0404             sqis=sq2/qr(4,iq);         <span class="comment">% inverse std dev of gaussian (proportional to slope)</span>
0405         <span class="keyword">end</span>
0406         <span class="keyword">switch</span> nr(10)   <span class="comment">% switch on model type</span>
0407             <span class="keyword">case</span> 1
0408                 <span class="keyword">for</span> i=1:length(ires)
0409                     j=ires(i);
0410                     r0=res(j,3)==0;
0411                     wq2=wq2+log(r0+(1-2*r0)*(guess+pscale*((1+exp(sqis*xq2'-xtstd-repmat(sqis,1,nxq)*res(j,2))).^(-1)))); <span class="comment">%  P(l | r,x)</span>
0412                 <span class="keyword">end</span>
0413             <span class="keyword">case</span> 2
0414                 <span class="keyword">for</span> i=1:length(ires)
0415                     j=ires(i);
0416                     r0=res(j,3)==0;
0417                     wq2=wq2+log(r0+(1-2*r0)*(guess+pscale*normcdf(repmat(sqis,1,nxq)*res(j,2)-sqis*xq2'+xtstd))); <span class="comment">%  P(l | r,x)</span>
0418                 <span class="keyword">end</span>
0419         <span class="keyword">end</span>
0420     <span class="keyword">else</span>                                                    <span class="comment">% possibly do quadratic interpolation</span>
0421         wq2=max(reshape(wqi,nsq,nxq)-max(wqi(:)),wfl);        <span class="comment">% turn into a normalized, clipped matrix for easy interpolation</span>
0422         <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0423         <span class="comment">% use quadratic interpolation in SRT axis %</span>
0424         <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0425         <span class="keyword">if</span> (xq2(end)-xq2(1))/(xqrange)&gt;nr(16)          <span class="comment">% if range has shrunk by &lt; nr(16), leave the SRT axis unchanged</span>
0426             xq2=xqi;                                  <span class="comment">% copy the old SRT axis</span>
0427             wqup=1;                                   <span class="comment">% update flag</span>
0428         <span class="keyword">else</span>                                                <span class="comment">% else do quadratic interpolation</span>
0429             <span class="keyword">if</span> LOG, fprintf(LOG,<span class="string">'N=%d:%d, interpolate x: x(%.2g,%.2g)-&gt;(%.2g,%.2g)\n'</span>,iq,nresq(iq),xqi([1 end]),xq2([1 end])); <span class="keyword">end</span>
0430             xqf=(xq2-xqi(1))/(xqi(2)-xqi(1));
0431             xqj=ceil(xqf);
0432             xqf=xqj-xqf;
0433             xqg=1-xqf;
0434             xqh=0.25*xqf.*xqg; <span class="comment">% Quadratic coefficient</span>
0435             xqf((xqj&lt;=0) | (xqj&gt;nxq))=0;
0436             xqg((xqj&lt;0) | (xqj&gt;=nxq))=0;
0437             xqh((xqj&lt;1) | (xqj&gt;=nxq-1))=0;
0438             wq2=wq2(:,min(max(xqj,1),nxq)).*repmat(xqf'+xqh',nsq,1)+wq2(:,min(max(xqj+1,1),nxq)).*repmat(xqg'+xqh',nsq,1) - <span class="keyword">...</span>
0439                 (wq2(:,min(max(xqj-1,1),nxq))+wq2(:,min(max(xqj+2,1),nxq))).*repmat(xqh',nsq,1);
0440         <span class="keyword">end</span>
0441         <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0442         <span class="comment">% use quadratic interpolation in slope axis %</span>
0443         <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0444         <span class="keyword">if</span> (sq2(end)-sq2(1))/(sqrange)&gt;nr(16)       <span class="comment">% if range has shrunk by &lt; nr(16), leave the slope axis unchanged</span>
0445             sq2=sqi;                                <span class="comment">% copy the old slope axis</span>
0446             wqup=wqup-1;                            <span class="comment">% update flag</span>
0447         <span class="keyword">else</span>
0448             <span class="keyword">if</span> LOG, fprintf(LOG,<span class="string">'N=%d:%d, interpolate s: s(%.2g,%.2g)-&gt;(%.2g,%.2g)\n'</span>,iq,nresq(iq),sqi([1 end]),sq2([1 end])); <span class="keyword">end</span>
0449             sqf=(sq2-sqi(1))/(sqi(2)-sqi(1));
0450             sqj=ceil(sqf);
0451             sqf=sqj-sqf;
0452             sqg=1-sqf;
0453             sqh=0.25*sqf.*sqg; <span class="comment">% Quadratic coefficient</span>
0454             sqf((sqj&lt;=0) | (sqj&gt;nsq))=0;
0455             sqg((sqj&lt;0) | (sqj&gt;=nsq))=0;
0456             sqh((sqj&lt;1) | (sqj&gt;=nsq-1))=0;
0457             wq2=wq2(min(max(sqj,1),nsq),:).*repmat(sqf+sqh,1,nxq)+wq2(min(max(sqj+1,1),nsq),:).*repmat(sqg+sqh,1,nxq) - <span class="keyword">...</span>
0458                 (wq2(min(max(sqj-1,1),nsq),:)+wq2(min(max(sqj+2,1),nsq),:)).*repmat(sqh,1,nxq);
0459         <span class="keyword">end</span>
0460     <span class="keyword">end</span>
0461     <span class="keyword">if</span> wqup&gt;0                               <span class="comment">% check if we have updated either axis</span>
0462         wq2=max(wq2(:)-max(wq2(:)),wfl);    <span class="comment">% turn back into a normalized, clipped vector</span>
0463         sqi=sq2;                            <span class="comment">% update slope (or log slope) values in PDF</span>
0464         xqi=xq2;                            <span class="comment">% update SRT values of the pdf array</span>
0465         wqi=wq2;                            <span class="comment">% log pdf array</span>
0466         sq(:,iq)=sqi;                       <span class="comment">% save new s axis (log slope or slope)</span>
0467         xq(:,iq)=xqi;                       <span class="comment">% save new x axis (SRT)</span>
0468         wq(:,iq)=wqi;                       <span class="comment">% save log-pdf</span>
0469     <span class="keyword">end</span>
0470     <span class="keyword">if</span> nr(11)                   <span class="comment">% if log slope</span>
0471         sqis=exp(sqi)/qr(4,iq);    <span class="comment">% inverse std dev of gaussian (proportional to slope)</span>
0472     <span class="keyword">else</span>
0473         sqis=sqi/qr(4,iq);         <span class="comment">% inverse std dev of gaussian (proportional to slope)</span>
0474     <span class="keyword">end</span>
0475     <span class="keyword">if</span> nargin&gt;=3
0476         <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0477         <span class="comment">% update pdfs with a new probe result  %</span>
0478         <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0479         <span class="keyword">if</span> nres&gt;=size(res,1)        <span class="comment">% ensure there is space for a new result</span>
0480             res=[res; zeros(nres,7)]; <span class="comment">% double the size of the results array</span>
0481         <span class="keyword">end</span>
0482         nres=nres+1;                <span class="comment">% increment the total number of results</span>
0483         nresq(iq)=nresq(iq)+1;        <span class="comment">% increment the number of results for this model</span>
0484         res(nres,1:3)=[iq,x,r];       <span class="comment">% save the latest result</span>
0485         xmm(:,iq)=[min(x,xmm(1,iq)); max(x,xmm(2,iq))]; <span class="comment">% update min/max probe values</span>
0486         hn=hn*hfact;         <span class="comment">% increase v_entropy gains to ensure all models get a chance</span>
0487         
0488         <span class="comment">%     disp(res(1:nres,:));  % ********** diagnostic</span>
0489         
0490         <span class="comment">%</span>
0491         <span class="comment">% update log probabilities with the previous test result</span>
0492         <span class="comment">% If r==1, we multiply (add in log) by a horizontally reflected version of the psychometric function that equals</span>
0493         <span class="comment">% the threshold (e.g. 0.5 or 0.75) at the probe SNR, x, for all slopes.</span>
0494         <span class="comment">% If r==0, we multiply (add in log) by 1- this reflected version which therefore equals (1-thresh) at x.</span>
0495         <span class="comment">%</span>
0496         r0=r==0;
0497         <span class="keyword">switch</span> nr(10)   <span class="comment">% switch on model type</span>
0498             <span class="keyword">case</span> 1
0499                 wqi=wqi+log(r0+(1-2*r0)*(guess+pscale*((1+exp(reshape(sqis*xqi'-xtstd,nsxq,1)-repmat(sqis,nxq,1)*x)).^(-1)))); <span class="comment">%  P(l | r,x)</span>
0500             <span class="keyword">case</span> 2
0501                 wqi=wqi+log(r0+(1-2*r0)*(guess+pscale*normcdf(repmat(sqis,nxq,1)*x-reshape(sqis*xqi'-xtstd,nsxq,1)))); <span class="comment">%  P(l | r,x)</span>
0502             <span class="keyword">otherwise</span>
0503                 error(<span class="string">'Unrecognised psychometric model selection'</span>);
0504         <span class="keyword">end</span>
0505         wq(:,iq)=wqi;             <span class="comment">% save updated probabilities</span>
0506     <span class="keyword">end</span>
0507     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0508     <span class="comment">% Calculate mean and covariance and v_entropy  %</span>
0509     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0510     ewqi=exp(wqi-max(wqi));             <span class="comment">% unnormalized probability vector</span>
0511     ewqi=ewqi/sum(ewqi);
0512     wqsx=reshape(ewqi,nsq,nxq);         <span class="comment">% normalized probabilities</span>
0513     px=sum(wqsx,1);                             <span class="comment">%  p(x0)</span>
0514     ps=sum(wqsx,2);                             <span class="comment">%  p(s0)</span>
0515     xe=px*xqi;                                  <span class="comment">% E(x0)</span>
0516     se=ps'*sqi;                                 <span class="comment">% E(s0)</span>
0517     
0518     [pxpk,xm]=max(px);                          <span class="comment">% marginal mode in x</span>
0519     <span class="keyword">if</span> xm&gt;1 &amp;&amp; xm&lt;nxq                           <span class="comment">% use quadratic interpolation in log prob if possible</span>
0520         [dum,xm2]=quadpeak(log(px(xm-1:xm+1))');
0521         xm=xm+xm2-2;
0522     <span class="keyword">end</span>
0523     xm=(2-xm)*xqi(1)+(xm-1)*xqi(2);             <span class="comment">% marginal mode(x0)</span>
0524     
0525     [pspk,sm]=max(ps);
0526     <span class="keyword">if</span> sm&gt;1 &amp;&amp; sm&lt;nsq                           <span class="comment">% use quadratic interpolation in log prob if possible</span>
0527         [dum,sm2]=<a href="v_quadpeak.html" class="code" title="function [v,x,t,m,ze]=v_quadpeak(z)">v_quadpeak</a>(log(ps(sm-1:sm+1)));
0528         sm=sm+sm2-2;
0529     <span class="keyword">end</span>
0530     sm=(2-sm)*sqi(1)+(sm-1)*sqi(2);             <span class="comment">% marginal mode(s0)</span>
0531     
0532     [wqpk,j]=max(ewqi);
0533     i=1+floor((j-1)/nsq);
0534     j=j-nsq*(i-1);
0535     <span class="keyword">if</span> i&gt;1 &amp;&amp; i&lt;nxq &amp;&amp; j&gt;1 &amp;&amp; j&lt;nsq             <span class="comment">% use quadratic interpolation in log prob if possible</span>
0536         [dum,ji]=quadpeak(wqi(repmat((j-1:j+1)',1,3)+repmat(nsq*(i-2:i),3,1)));
0537         i=i+ji(2)-2;
0538         j=j+ji(1)-2;
0539     <span class="keyword">end</span>
0540     xj=(2-i)*xqi(1)+(i-1)*xqi(2);               <span class="comment">% joint mode  x</span>
0541     sj=(2-j)*sqi(1)+(j-1)*sqi(2);               <span class="comment">% joint mode: s</span>
0542     xv=px*(xqi.^2)-xe^2;                        <span class="comment">% Var(x0)</span>
0543     sv=ps'*(sqi.^2)-se^2;                       <span class="comment">% Var(s0)</span>
0544     sxv=ewqi'*(repmat(sqi-se,nxq,1).*reshape(repmat(xqi'-xe,nsq,1),nsxq,1)); <span class="comment">% Cov(s0*x0)</span>
0545     mq(:,iq,1)=[xe; se];                        <span class="comment">% save means</span>
0546     mq(:,iq,2)=[xj; sj];                        <span class="comment">% save joint mode</span>
0547     mq(:,iq,3)=[xm; sm];                        <span class="comment">% save marginal modes</span>
0548     vq(:,iq)=[xv; sxv; sv];                     <span class="comment">% save covariance matrix</span>
0549     xh=(px*log(px)')*(xqi(1)-xqi(2));           <span class="comment">% differential v_entropy h(x0)</span>
0550     sh=(ps'*log(ps))*(sqi(1)-sqi(2));              <span class="comment">% differential v_entropy h(s0)</span>
0551     <span class="keyword">if</span> nargin==3
0552         <span class="keyword">switch</span> nr(9) <span class="comment">% cost function: 1=variance, 2=v_entropy</span>
0553             <span class="keyword">case</span> 1
0554                 res(nres,4:7)=[xe se xv sv];               <span class="comment">% save info for plotting history</span>
0555             <span class="keyword">case</span> 2
0556                 res(nres,4:7)=[xe se xh sh];                                           <span class="comment">% find the minimum of cost function</span>
0557         <span class="keyword">end</span>
0558     <span class="keyword">end</span>
0559     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0560     <span class="comment">% now calculate robust mean and covariance for all models   %</span>
0561     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0562     <span class="keyword">if</span> nargout&gt;4 || ~isempty(po)       <span class="comment">% need to update robust statistics if &gt;4 outputs or if plotting</span>
0563         <span class="keyword">if</span>  nr(19)==0               <span class="comment">% if no outliers allowed in any case</span>
0564             mr=mq;
0565             vr=vq;
0566         <span class="keyword">elseif</span> nresr==nres          <span class="comment">% if robust statistics are all up-to-date</span>
0567             mr=mqr;
0568             vr=vqr;
0569         <span class="keyword">else</span>                        <span class="comment">% update robust statistics for all models</span>
0570             <span class="keyword">for</span> jq=1:length(nresq)             <span class="comment">% loop for each model</span>
0571                 <span class="keyword">if</span> any(res(nresr+1:nres,1)==jq)   <span class="comment">% if we have had any recent probe results for this model</span>
0572                     guessj=pr(3,jq);                <span class="comment">% guess rate (1/choices)</span>
0573                     pscalej=qr(1,jq);            <span class="comment">% prob range left after subtracting miss and guess probs</span>
0574                     xtstdj=qr(3,jq);               <span class="comment">% x position of target in std measure</span>
0575                     xqj=xq(:,jq);               <span class="comment">% SRT values of the pdf array</span>
0576                     sqj=sq(:,jq);               <span class="comment">% slope (or log slope) values in PDF</span>
0577                     <span class="comment">% first determine the SNR range to include</span>
0578                     xej=mq(1,jq,1);  <span class="comment">% mean SRT</span>
0579                     sej=mq(2,jq,1);   <span class="comment">% mean slope or log slope</span>
0580                     resj=res(res(:,1)==jq,:);  <span class="comment">% results for this model</span>
0581                     <span class="keyword">if</span> nr(11)                   <span class="comment">% if log slope</span>
0582                         sqisj=exp(sej)/qr(4,jq);    <span class="comment">% inverse std dev of gaussian (proportional to slope) at mean slope</span>
0583                     <span class="keyword">else</span>
0584                         sqisj=sej/qr(4,jq);         <span class="comment">% inverse std dev of gaussian (proportional to slope) at mean slope</span>
0585                     <span class="keyword">end</span>
0586                     <span class="comment">% find upper bound for SNR with r=0</span>
0587                     <span class="keyword">switch</span> nr(10)   <span class="comment">% switch on model type</span>
0588                         <span class="keyword">case</span> 1
0589                             xlim(2,jq)=xej-(xtstdj+log(nr(19)/(pscalej-nr(19))))/sqisj; <span class="comment">% outlier if r=0 and x &gt; xlim</span>
0590                         <span class="keyword">case</span> 2
0591                             xlim(2,jq)=xej+(norminv((pscalej-nr(19))/pscalej)-xtstdj)/sqisj; <span class="comment">% outlier if r=1 and x &lt; xlim</span>
0592                     <span class="keyword">end</span>
0593                     mou=resj(:,2)&gt;=xlim(2,jq) &amp; resj(:,3)==0;       <span class="comment">% find high outliers</span>
0594                     <span class="keyword">if</span> nr(19)&gt;guessj  <span class="comment">% find lower bound for SNR with r=1</span>
0595                         <span class="keyword">switch</span> nr(10)   <span class="comment">% switch on model type</span>
0596                             <span class="keyword">case</span> 1
0597                                 xlim(1,jq)=xej-(xtstdj+log(pscalej/(nr(19)-guessj)-1))/sqisj; <span class="comment">% outlier if r=1 and x &lt; xlim</span>
0598                             <span class="keyword">case</span> 2
0599                                 xlim(1,jq)=xej+(norminv((nr(19)-guessj)/pscalej)-xtstdj)/sqisj; <span class="comment">% outlier if r=1 and x &lt; xlim</span>
0600                         <span class="keyword">end</span>
0601                         mou=mou | resj(:,2)&lt;=xlim(1,jq) &amp; resj(:,3)==1;       <span class="comment">% add in low outliers</span>
0602                     <span class="keyword">end</span>
0603                     <span class="keyword">if</span> any(mou) <span class="comment">% if any results must be excluded</span>
0604                         <span class="comment">% now recalculate pdf excluding bad points</span>
0605                         <span class="keyword">if</span> nr(11)                   <span class="comment">% if log slope</span>
0606                             sqisj=exp(sqj)/qr(4,jq);    <span class="comment">% inverse std dev of gaussian (proportional to slope) for all grid points</span>
0607                         <span class="keyword">else</span>
0608                             sqisj=sqj/qr(4,jq);         <span class="comment">% inverse std dev of gaussian (proportional to slope) for all grid points</span>
0609                         <span class="keyword">end</span>
0610                         wqj=repmat(pq0(2,jq)*(sqj-mq0(2,jq)).^2,1,nxq)+repmat(pq0(1,jq)*(xqj-mq0(1,jq))'.^2,nsq,1); <span class="comment">% log probs of prior</span>
0611                         <span class="keyword">switch</span> nr(10)               <span class="comment">% switch on model type</span>
0612                             <span class="keyword">case</span> 1                  <span class="comment">% nr(10)=1: logistic model</span>
0613                                 <span class="keyword">for</span> j=find(mou==0)'
0614                                     r0=resj(j,3)==0;
0615                                     wqj=wqj+log(r0+(1-2*r0)*(guessj+pscalej*((1+exp(sqisj*xqj'-xtstdj-repmat(sqisj,1,nxq)*resj(j,2))).^(-1)))); <span class="comment">%  P(l | r,x)</span>
0616                                 <span class="keyword">end</span>
0617                             <span class="keyword">case</span> 2                  <span class="comment">% nr(10)=2: cumulative gaussian model</span>
0618                                 <span class="keyword">for</span> j=find(mou==0)'
0619                                     r0=resj(j,3)==0;
0620                                     wqj=wqj+log(r0+(1-2*r0)*(guessj+pscalej*normcdf(repmat(sqisj,1,nxq)*resj(j,2)-sqisj*xqj'+xtstdj))); <span class="comment">%  P(l | r,x)</span>
0621                                 <span class="keyword">end</span>
0622                         <span class="keyword">end</span>
0623                         <span class="comment">% now calculate robust means and modes</span>
0624                         wqj=wqj(:);                         <span class="comment">% turn 2D array into a column vector</span>
0625                         ewqj=exp(wqj-max(wqj));             <span class="comment">% unnormalized probability vector</span>
0626                         ewqj=ewqj/sum(ewqj);
0627                         wqsxr=reshape(ewqj,nsq,nxq);         <span class="comment">% normalized probabilities</span>
0628                         pxr=sum(wqsxr,1);                             <span class="comment">%  p(x0)</span>
0629                         psr=sum(wqsxr,2);                             <span class="comment">%  p(s0)</span>
0630                         xer=pxr*xqj;                                  <span class="comment">% E(x0)</span>
0631                         ser=psr'*sqj;                                 <span class="comment">% E(s0)</span>
0632                         [pxpk,xmr]=max(pxr);                          <span class="comment">% marginal mode in x</span>
0633                         <span class="keyword">if</span> xmr&gt;1 &amp;&amp; xmr&lt;nxq                           <span class="comment">% use quadratic interpolation in log prob if possible</span>
0634                             [dum,xm2]=quadpeak(log(pxr(xmr-1:xmr+1))');
0635                             xmr=xmr+xm2-2;
0636                         <span class="keyword">end</span>
0637                         xmr=(2-xmr)*xqj(1)+(xmr-1)*xqj(2);             <span class="comment">% marginal mode(x0) (requires the xqi to be uniformly spaced)</span>
0638                         [pspk,smr]=max(psr);
0639                         <span class="keyword">if</span> smr&gt;1 &amp;&amp; smr&lt;nsq                           <span class="comment">% use quadratic interpolation in log prob if possible</span>
0640                             [dum,sm2]=<a href="v_quadpeak.html" class="code" title="function [v,x,t,m,ze]=v_quadpeak(z)">v_quadpeak</a>(log(psr(smr-1:smr+1)));
0641                             smr=smr+sm2-2;
0642                         <span class="keyword">end</span>
0643                         smr=(2-smr)*sqj(1)+(smr-1)*sqj(2);            <span class="comment">% marginal mode(s0)</span>
0644                         [wqpk,j]=max(ewqj);                         <span class="comment">% find max of 1-D array</span>
0645                         i=1+floor((j-1)/nsq);                       <span class="comment">% convert to 2-D index (i,j)</span>
0646                         j=j-nsq*(i-1);
0647                         <span class="keyword">if</span> i&gt;1 &amp;&amp; i&lt;nxq &amp;&amp; j&gt;1 &amp;&amp; j&lt;nsq             <span class="comment">% use quadratic interpolation in log prob if possible</span>
0648                             [dum,ji]=quadpeak(wqj(repmat((j-1:j+1)',1,3)+repmat(nsq*(i-2:i),3,1)));
0649                             i=i+ji(2)-2;
0650                             j=j+ji(1)-2;
0651                         <span class="keyword">end</span>
0652                         xjr=(2-i)*xqj(1)+(i-1)*xqj(2);               <span class="comment">% joint mode  x</span>
0653                         sjr=(2-j)*sqj(1)+(j-1)*sqj(2);               <span class="comment">% joint mode: s</span>
0654                         xvr=pxr*(xqj.^2)-xer^2;                       <span class="comment">% Var(x0)</span>
0655                         svr=psr'*(sqj.^2)-ser^2;                      <span class="comment">% Var(s0)</span>
0656                         sxvr=ewqj'*(repmat(sqj-ser,nxq,1).*reshape(repmat(xqj'-xer,nsq,1),nsxq,1)); <span class="comment">% Cov(s0*x0)</span>
0657                         mqr(:,jq,:)=reshape([xer xjr xmr; ser sjr smr],[2 1 3]); <span class="comment">% save means, joint modes and marginal modes</span>
0658                         vqr(:,jq)=[xvr; sxvr; svr];                    <span class="comment">% save covariance matrix</span>
0659                     <span class="keyword">else</span>   <span class="comment">% if there are no outliers, just copy non-robust results</span>
0660                         mqr(:,jq,:)= mq(:,jq,:);                    <span class="comment">% use non-robust means</span>
0661                         vqr(:,jq)= vq(:,jq);                          <span class="comment">% and non-robust variances</span>
0662                     <span class="keyword">end</span>
0663                 <span class="keyword">end</span>
0664             <span class="keyword">end</span>
0665             nresr=nres;                 <span class="comment">% all robust means/modes are now up-to-date</span>
0666             mr=mqr;                     <span class="comment">% send to output arguments</span>
0667             vr=vqr;
0668         <span class="keyword">end</span>
0669         <span class="keyword">if</span> nr(11)    <span class="comment">% if using log-slope</span>
0670             mr(2,:,:)=exp(mr(2,:,:));  <span class="comment">% convert to real slope</span>
0671             vr(2,:)=vr(2,:).*mr(2,:,1); <span class="comment">% correct the covariance</span>
0672             vr(3,:)=vr(3,:).*mr(2,:,1).^2; <span class="comment">% and the slope variance</span>
0673         <span class="keyword">end</span>
0674     <span class="keyword">end</span>
0675     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0676     <span class="comment">% now estimate the next probe SNR            %</span>
0677     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0678     <span class="keyword">if</span> ~numel(xz{iq})                               <span class="comment">% if no list of probe SNRs was specified</span>
0679         ytry=exp(0.25i*pi*(0:7))';                  <span class="comment">% points around the circle</span>
0680         ytry=[real(ytry) imag(ytry)];
0681         [vtry,dtry]=eig([xv sxv; sxv sv]);          <span class="comment">% eigendecomposition of covariance matrix</span>
0682         tryxs=repmat([xe,se],8,1)+nr(14)*ytry*sqrt(dtry)*vtry';
0683         pmin=0.05;                                  <span class="comment">% target probe success probability</span>
0684         <span class="keyword">if</span> nr(11)
0685             tryxs(:,2)=qr(4,iq)*exp(-tryxs(:,2));     <span class="comment">% convert log(slope) to std dev</span>
0686         <span class="keyword">else</span>
0687             tryxs(:,2)=qr(4,iq)*tryxs(:,2).^(-1);     <span class="comment">% convert slope to std dev</span>
0688         <span class="keyword">end</span>
0689         dp=nr(15);  <span class="comment">% maximum shift of probe value outside previous range</span>
0690         <span class="keyword">switch</span> nr(10)   <span class="comment">% switch by model type</span>
0691             <span class="keyword">case</span> 1 <span class="comment">% logistic</span>
0692                 qmax=min(xmm(2,iq)+dp,max(tryxs(:,1)+(log((1-pmin)/pmin)-xtstd)*tryxs(:,2)));
0693                 qmin=max(xmm(1,iq)-dp,min(tryxs(:,1)+(log(pmin/(1-pmin))-xtstd)*tryxs(:,2)));
0694             <span class="keyword">case</span> 2 <span class="comment">% cumulative gaussian</span>
0695                 qmax=min(xmm(2,iq)+dp, max(tryxs(:,1)+(norminv(1-pmin)-xtstd)*tryxs(:,2)));
0696                 qmin=max(xmm(1,iq)-dp,min(tryxs(:,1)+(norminv(pmin)-xtstd)*tryxs(:,2)));
0697         <span class="keyword">end</span>
0698         dxt=max(nr(5),(qmax-qmin)/nxh);             <span class="comment">% minimum step size of nr(5) [0.2 dB]</span>
0699         xt=(qmin+qmax)/2+((1:nxh)-(1+nxh)/2)*dxt;
0700     <span class="keyword">else</span>                                            <span class="comment">% if a specific list of probe SNRs exists</span>
0701         xzi=xz{iq};                                 <span class="comment">% xzi is the list of available probe SNRs</span>
0702         <span class="keyword">if</span> numel(xzi)&lt;=nxh                          <span class="comment">% use all available probe SNRs if there are not too many</span>
0703             xt=xzi;
0704         <span class="keyword">else</span>
0705             [xt,ixt]=min(abs(xzi-xe));                                  <span class="comment">% find the closest one to xe ** not necessarily optimum ***</span>
0706             ixt=max(1,min(1+numel(xzi)-nxh,ixt-floor((1+nxh)/2)));      <span class="comment">% arrange symmetrically around xt</span>
0707             xt=xzi(ixt:min(ixt+nxh-1,numel(xzi)));
0708         <span class="keyword">end</span>
0709     <span class="keyword">end</span>
0710     nxhp=length(xt);  <span class="comment">% xt are the potential probe SNRs</span>
0711     <span class="comment">% Now find the probe value that minimizes the expected value of the cost function</span>
0712     <span class="comment">% In the following: l = parameters of psychometric function, x = the probe SNR and r = the probe result</span>
0713     <span class="keyword">switch</span> nr(10)
0714         <span class="keyword">case</span> 1
0715             prt=guess+pscale*((1+exp(repmat(reshape(sqis*xqi'-xtstd,nsxq,1),1,nxhp)-repmat(sqis,nxq,1)*xt)).^(-1)); <span class="comment">%  P(r=1 | l,x)</span>
0716         <span class="keyword">case</span> 2
0717             prt=guess+pscale*normcdf(repmat(sqis,nxq,1)*xt-repmat(reshape(sqis*xqi'-xtstd,nsxq,1),1,nxhp)); <span class="comment">%  P(r=1 | l,x)</span>
0718     <span class="keyword">end</span>
0719     wqt=repmat(ewqi,1,nxhp);
0720     hminr=zeros(2,1);                   <span class="comment">% space for minimum expected cost function for each r0</span>
0721     hminj=zeros(1,nxhp);   <span class="comment">% space for minimum expected cost function for each x0</span>
0722     <span class="keyword">if</span> nr(18) <span class="comment">% if doing look 2-ahead</span>
0723         pl1=prt.*wqt;                       <span class="comment">% posterior prob given success = p(l | x0,r0=1) unnormalized</span>
0724         pl0=wqt-pl1;                        <span class="comment">% posterior prob given failure = p(l | x0,r0=0) unnormalized</span>
0725         pr0=sum(pl1,1);                     <span class="comment">% p(r0=1 | x0)=Sum{P(r0=1 | l,x0)*P(l)} [note each column of wqt is normalized] (row vector)</span>
0726         plxr=reshape([pl0./repmat(1-pr0,nsxq,1) pl1./repmat(pr0,nsxq,1)],nsxq,nxhp,2); <span class="comment">% posterior prob p(l | x0,r0) column-normalized</span>
0727         nx0=nxhp; <span class="comment">% outer loop for each x0</span>
0728         nr0=2; <span class="comment">% inner loop for each r0</span>
0729     <span class="keyword">else</span>  <span class="comment">% if only doing look 1-ahead</span>
0730         nx0=1;  <span class="comment">% only execute outer loop once</span>
0731         nr0=1;  <span class="comment">% only execute inner loop once</span>
0732         pr0=0; <span class="comment">% dummy value (used at end of outer loop)</span>
0733     <span class="keyword">end</span>
0734     hx2=zeros(nx0,nxhp,nr0);        <span class="comment">% space for square array of expected cost functions (for possible plotting only)</span>
0735     <span class="keyword">for</span> j=1:nx0                     <span class="comment">% loop for each possible probe SNR, x0, (or only once is look-1-ahead)</span>
0736         <span class="keyword">for</span> jr=1:nr0                 <span class="comment">% loop for each possible test result, r0=jr-1 (or only once is look-1-ahead)</span>
0737             <span class="keyword">if</span> nr(18) <span class="comment">% if doing look 2-ahead</span>
0738                 wqt=repmat(plxr(:,j,jr),1,nxhp);         <span class="comment">% posterior prob p(l | x0=xt(j),r0=jr-1) column-normalized</span>
0739             <span class="keyword">end</span>
0740             pl1=prt.*wqt;                     <span class="comment">% posterior prob given success = p(l | x,r=1) unnormalized</span>
0741             pl0=wqt-pl1;                      <span class="comment">% posterior prob given failure = p(l | x,r=0) unnormalized</span>
0742             prh=sum(pl1,1);                     <span class="comment">% p(r | x)=Sum{P(r | l,x)*P(l)} [note each column of wqtjr is normalized] (row vector)</span>
0743             pl1=pl1./repmat(prh,nsxq,1);        <span class="comment">% posterior prob given success = p(l | x,r=1) normalized</span>
0744             pl0=pl0./repmat(1-prh,nsxq,1);      <span class="comment">% posterior prob given failure = p(l | x,r=0) normalized</span>
0745             
0746             px1=squeeze(sum(reshape(pl1,nsq,nxq,[]),1));    <span class="comment">% p(x0 | x,r=1)</span>
0747             px0=squeeze(sum(reshape(pl0,nsq,nxq,[]),1));    <span class="comment">% p(x0 | x,r=0)</span>
0748             ps1=squeeze(sum(reshape(pl1,nsq,nxq,[]),2));    <span class="comment">% p(s0 | x,r=1)</span>
0749             ps0=squeeze(sum(reshape(pl0,nsq,nxq,[]),2));    <span class="comment">% p(s0 | x,r=0)</span>
0750             xet1=xqi'*px1;                                  <span class="comment">% E(x0 | x,r=1)</span>
0751             xvt1=(xqi.^2)'*px1-xet1.^2;                     <span class="comment">% Var(x0 | x,r=1)</span>
0752             xet0=xqi'*px0;                                  <span class="comment">% E(x0 | x,r=0)</span>
0753             xvt0=(xqi.^2)'*px0-xet0.^2;                     <span class="comment">% Var(x0 | x,r=0)</span>
0754             xvt=xvt1.*prh+xvt0.*(1-prh);                    <span class="comment">% E(Var(x0 | x ))</span>
0755             set1=sqi'*ps1;                                  <span class="comment">% E(s0 | x,r=1)</span>
0756             svt1=(sqi.^2)'*ps1-set1.^2;                     <span class="comment">% Var(s0 | x,r=1)</span>
0757             set0=sqi'*ps0;                                  <span class="comment">% E(s0 | x,r=0)</span>
0758             svt0=(sqi.^2)'*ps0-set0.^2;                     <span class="comment">% Var(s0 | x,r=0)</span>
0759             svt=svt1.*prh+svt0.*(1-prh);                    <span class="comment">% E(Var(s0 | x ))</span>
0760             xht1=sum(log(px1).*px1,1);                      <span class="comment">% -H(x0 | x,r=1)</span>
0761             xht0=sum(log(px0).*px0,1);                      <span class="comment">% -H(x0 | x,r=0)</span>
0762             xht=(xht1.*prh+xht0.*(1-prh))*(xqi(1)-xqi(2));    <span class="comment">% h(x0 | x)</span>
0763             sht1=sum(log(ps1).*ps1,1);                      <span class="comment">% -H(s0 | x,r=1)</span>
0764             sht0=sum(log(ps0).*ps0,1);                      <span class="comment">% -H(s0 | x,r=0)</span>
0765             sht=(sht1.*prh+sht0.*(1-prh))*(sqi(1)-sqi(2));    <span class="comment">% h(s0 | x)</span>
0766             <span class="keyword">switch</span> nr(9) <span class="comment">% cost function: 1=variance, 2=v_entropy</span>
0767                 <span class="keyword">case</span> 1
0768                     hx=(xvt + nr(4)*svt)/(1+nr(4));               <span class="comment">% expected cost function for each possible test SNR</span>
0769                 <span class="keyword">case</span> 2
0770                     hx=(xht + nr(4)*sht)/(1+nr(4));                <span class="comment">% expected cost function for each possible test SNR                                          % find the minimum of cost function</span>
0771             <span class="keyword">end</span>
0772             hx2(j,:,jr)=hx; <span class="comment">% save expected cost function for possible plotting</span>
0773             [hminr(jr),ix]=min(hx);       <span class="comment">% find the minimum of cost function for each value of r (0 or 1)</span>
0774             <span class="comment">%         fprintf('Probe range: %.3g %.3g; choose %.3g\n',xt(1),xt(end),xt(ix));</span>
0775         <span class="keyword">end</span>
0776         hminj(j)=[1-pr0(j) pr0(j)]*hminr; <span class="comment">% expected cost function for each x0 assuming x1 is chosen optimally</span>
0777     <span class="keyword">end</span>
0778     <span class="keyword">if</span> nr(18) <span class="comment">% if doing look 2-ahead</span>
0779         [hminr(1),ix]=min(hminj);                         <span class="comment">% find the minimum of cost function</span>
0780     <span class="keyword">end</span>
0781     <span class="comment">% now ix indexes the optimal probe choice and hminr(1) is the resultant expected cost function</span>
0782     xn(iq)=xt(ix);  <span class="comment">% save the next probe snr for this model</span>
0783     <span class="comment">% calculate the expected decrease in cost function (to decide which model to probe next)</span>
0784     <span class="keyword">switch</span> nr(9)  <span class="comment">% cost function: 1=variance, 2=v_entropy</span>
0785         <span class="keyword">case</span> 1                                              <span class="comment">% minimize variance</span>
0786             hn(iq)=(xv + nr(4)*sv)/(1+nr(4))-hminr(1);        <span class="comment">% save the expected decrease in cost function for this model</span>
0787         <span class="keyword">case</span> 2                                              <span class="comment">% minimize v_entropy</span>
0788             hn(iq)=(xh + nr(4)*sh)/(1+nr(4))-hminr(1);        <span class="comment">% save the expected decrease in cost function for this model</span>
0789     <span class="keyword">end</span>
0790     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0791     <span class="comment">%  now do plotting                                 %</span>
0792     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0793     <span class="keyword">if</span> ~isempty(po)
0794         h=gcf;
0795         <span class="keyword">if</span> isreal(h)
0796             fig0=round(h); <span class="comment">% get figure number</span>
0797         <span class="keyword">else</span>
0798             fig0=get(h,<span class="string">'number'</span>);  <span class="comment">% in new versions of matlab use this method</span>
0799         <span class="keyword">end</span>
0800         axm=[3 -2; 2 -1]; <span class="comment">% matrix to extend an axis beackwards by two steps</span>
0801         <span class="keyword">for</span> i=1:length(po)
0802             poi=po(i);
0803             <span class="keyword">switch</span> poi
0804                 <span class="keyword">case</span> <span class="string">'p'</span>    <span class="comment">% plot pdf</span>
0805                     nxq=nr(1);
0806                     nsq=nr(2);
0807                     xqj=[axm*xqi(1:2);xqi];
0808                     sqj=[axm*sqi(1:2);sqi];
0809                     figure(fig0);
0810                     imagesc(xqj,sqj,[zeros(1,2) px/pxpk; zeros(1,nxq+2);ps/pspk zeros(nsq,1) wqsx/wqpk]);
0811                     hold on
0812                     plot(xe,se,<span class="string">'ko'</span>,xm,sm,<span class="string">'k^'</span>,xj,sj,<span class="string">'k*'</span>);
0813                     plot(xqj(2),se,<span class="string">'wo'</span>,xqj(2),sm,<span class="string">'w^'</span>,xqj(2),sj,<span class="string">'w*'</span>);
0814                     plot(xe,sqj(2),<span class="string">'wo'</span>,xm,sqj(2),<span class="string">'w^'</span>,xj,sqj(2),<span class="string">'w*'</span>);
0815                     <span class="keyword">if</span> any(mqr(:,iq,:)~=mq(:,iq,:))
0816                         plot(mqr(1,iq,1),mqr(2,iq,1),<span class="string">'go'</span>,mqr(1,iq,3),mqr(2,iq,3),<span class="string">'g^'</span>,xj,sj,<span class="string">'g*'</span>);
0817                         plot(xqj(end),mqr(2,iq,1),<span class="string">'go'</span>,xqj(end),mqr(2,iq,3),<span class="string">'g^'</span>,xqj(end),mqr(2,iq,2),<span class="string">'g*'</span>);
0818                         plot(mqr(1,iq,1),sqj(end),<span class="string">'go'</span>,mqr(1,iq,3),sqj(end),<span class="string">'g^'</span>,mqr(1,iq,2),sqj(end),<span class="string">'g*'</span>);
0819                         title(<span class="string">'Joint pdf: o mean, * mode, \Delta marg mode (green=robust)'</span>);
0820                     <span class="keyword">else</span>
0821                         title(<span class="string">'Joint pdf: o mean, * mode, \Delta marginal mode'</span>);
0822                     <span class="keyword">end</span>
0823                     t=linspace(0,2*pi,200);
0824                     xcir=cos(t);
0825                     scir=sin(t);
0826                     vcir=sqrt((sv*xcir.^2+xv*scir.^2-2*sxv*xcir.*scir)/(xv*sv-sxv^2));
0827                     plot(xe+xcir./vcir,se+scir./vcir,<span class="string">'w-'</span>);
0828                     hold off
0829                     axis <span class="string">'xy'</span>;
0830                     <span class="comment">%     colorbar;</span>
0831                     <span class="comment">%     v_cblabel('Relative Probability');</span>
0832                     xlabel(sprintf(<span class="string">'SNR @ %d%%SRT (dB)'</span>,round(pr(1)*100)));
0833                     <span class="keyword">if</span> nr(11)
0834                         ylabel(<span class="string">'Log psychometric Slope at threshold (prob/dB)'</span>);
0835                     <span class="keyword">else</span>
0836                         ylabel(<span class="string">'Psychometric Slope at threshold (prob/dB)'</span>);
0837                     <span class="keyword">end</span>
0838                     fig0=fig0+1;
0839                 <span class="keyword">case</span> <span class="string">'c'</span> <span class="comment">% plot cost evolution</span>
0840                     resi=res(res(:,1)==iq,:);
0841                     nresi=size(resi,1);
0842                     <span class="keyword">if</span> nresi&gt;0
0843                         figure(fig0);
0844                         plot(1:nresi,resi(:,6)/(1+nr(4)),<span class="string">'--r'</span>,1:nresi,(resi(:,6) + nr(4)*resi(:,6))/(1+nr(4)),<span class="string">'-b'</span>);
0845                         <a href="v_axisenlarge.html" class="code" title="function v_axisenlarge(f,h)">v_axisenlarge</a>([-1 -1.03]);
0846                         xlabel(<span class="string">'Trial Number'</span>);
0847                         <span class="keyword">switch</span> nr(9)  <span class="comment">% cost function: 1=variance, 2=v_entropy</span>
0848                             <span class="keyword">case</span> 1                                              <span class="comment">% minimize variance</span>
0849                                 ylabel(<span class="string">'Weighted Variance: SRT and Total'</span>);
0850                             <span class="keyword">case</span> 2                                              <span class="comment">% minimize v_entropy</span>
0851                                 ylabel(<span class="string">'Weighted Entropy: SRT and Total'</span>);
0852                         <span class="keyword">end</span>
0853                         title (<span class="string">'Cost function evolution'</span>);
0854                         fig0=fig0+1;
0855                     <span class="keyword">end</span>
0856                 <span class="keyword">case</span> <span class="string">'h'</span> <span class="comment">% plot history</span>
0857                     resi=res(res(:,1)==iq,:);
0858                     nresi=size(resi,1);
0859                     <span class="keyword">if</span> nresi&gt;0
0860                         <span class="keyword">if</span> nr(11)
0861                             resi(:,5)=exp(resi(:,5)); <span class="comment">% convert to real slope</span>
0862                             ylim=xe+([-2 0 thresh 1 3]-thresh)/exp(se);
0863                         <span class="keyword">else</span>
0864                             ylim=xe+([-2 0 thresh 1 3]-thresh)/se;
0865                         <span class="keyword">end</span>
0866                         figure(fig0);
0867                         plot(1:nresi,resi(:,[4 4 4])+resi(:,5).^(-1)*([0 thresh 1]-thresh),<span class="string">'-b'</span>);
0868                         hold on;
0869                         <span class="comment">% could plot outliers in a different colour by checking mou(res(:,1)==iq)</span>
0870                         msky=resi(:,3)==1;
0871                         plot(find(msky),resi(msky,2),<span class="string">'r+'</span>);
0872                         msky=resi(:,3)==0;
0873                         plot(find(msky),resi(msky,2),<span class="string">'bo'</span>);
0874                         gcaylim=get(gca,<span class="string">'ylim'</span>);
0875                         poutb=0;
0876                         <span class="keyword">if</span> nr(19)&gt;0 &amp;&amp; xlim(2,iq)&lt;gcaylim(2) <span class="comment">% plot upper outlier bound</span>
0877                             plot([1;nresi],[1;1]*xlim(2,iq),<span class="string">'-r'</span>);
0878                             <a href="v_texthvc.html" class="code" title="function h=v_texthvc(x,y,t,p,q,r)">v_texthvc</a>(nresi,xlim(2,iq),<span class="string">'max o'</span>,<span class="string">'rtr'</span>);
0879                         <span class="keyword">end</span>
0880                         <span class="keyword">if</span> nr(19)&gt;guess &amp;&amp; xlim(1,iq)&gt;gcaylim(1) <span class="comment">% plot upper outlier bound</span>
0881                             plot([1;nresi],[1;1]*xlim(1,iq),<span class="string">'-r'</span>);
0882                             <a href="v_texthvc.html" class="code" title="function h=v_texthvc(x,y,t,p,q,r)">v_texthvc</a>(nresi,xlim(1,iq),<span class="string">'min+'</span>,<span class="string">'rbr'</span>);
0883                         <span class="keyword">end</span>
0884                         hold off;
0885                         xlabel(<span class="string">'Trial Number'</span>);
0886                         ylabel(<span class="string">'Probe SNR (dB) [o,x=result]'</span>);
0887                         <a href="v_axisenlarge.html" class="code" title="function v_axisenlarge(f,h)">v_axisenlarge</a>(-1.02);
0888                         title(<span class="string">'SRT mean + slope intercepts'</span>);
0889                         fig0=fig0+1;
0890                     <span class="keyword">end</span>
0891                 <span class="keyword">case</span> <span class="string">'x'</span> <span class="comment">% plot probe choice</span>
0892                     figure(fig0);
0893                     <span class="keyword">if</span> nr(18) <span class="comment">% if doing look 2-ahead</span>
0894                         xtaxm=xt(1:2)*axm';
0895                         hx2e=hx2(:,:,1).*repmat(1-pr0',1,nxhp)+hx2(:,:,2).*repmat(pr0',1,nxhp);
0896                         [dum,ix2m]=min(hx2,[],2); <span class="comment">% ix2 indices to minimize hx2</span>
0897                         <span class="comment">% ==========  note that the following line assumes that the xt()</span>
0898                         <span class="comment">% are uniformly spaced (not necessarily true)</span>
0899                         imagesc([xtaxm xt],xt,[hminj' repmat(max(hx2e(:)),nxhp,1) hx2e]);
0900                         axis xy;
0901                         colorbar;
0902                         hold on
0903                         [hnmin,ii]=max(hn);         <span class="comment">% find the next model to probe</span>
0904                         plot(xtaxm(1),xn(iq),<span class="string">'&gt;w'</span>);
0905                         <span class="keyword">if</span> iq==ii
0906                             plot(xt(1:2)*axm(1,:)',xn(iq),<span class="string">'*w'</span>); <span class="comment">% circle the next model to be probed</span>
0907                         <span class="keyword">end</span>
0908                         plot([xtaxm(1) xt(end)],xn(iq)*[1 1],<span class="string">':w'</span>);
0909                         plot(xt(ix2m(:,1,1)),xt,<span class="string">'ow'</span>);
0910                         plot(xt(ix2m(:,1,2)),xt,<span class="string">'+w'</span>);
0911                         hold off
0912                         xlabel(<span class="string">'SNR probe 2 (dB) [+,o=best for R1=1,0]'</span>);
0913                         ylabel(sprintf(<span class="string">'SNR probe 1 (dB) [&gt; = %.2f]'</span>,xn(iq)));
0914                         title(<span class="string">'Expected Cost after 2 probes'</span>);
0915                     <span class="keyword">else</span> <span class="comment">% if doing look 1-ahead only</span>
0916                         plot(xt,hx2(1,:,1),<span class="string">'b'</span>);
0917                         <a href="v_axisenlarge.html" class="code" title="function v_axisenlarge(f,h)">v_axisenlarge</a>([-1 -1.03]);
0918                         caylim=get(gca,<span class="string">'ylim'</span>);
0919                         hold on
0920                         plot(xn([iq iq]),caylim,<span class="string">':k'</span>);
0921                         xlabel(sprintf(<span class="string">'Next probe SNR (dB) [min @ %.2f]'</span>,xn(iq)));
0922                         ylabel(<span class="string">'Cost'</span>)
0923                         title(<span class="string">'Expected Cost after next probe'</span>);
0924                     <span class="keyword">end</span>
0925                     fig0=fig0+1;
0926             <span class="keyword">end</span>
0927         <span class="keyword">end</span>
0928     <span class="keyword">end</span>
0929 <span class="keyword">end</span>
0930 <span class="keyword">if</span> iq~=0
0931     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0932     <span class="comment">%  now select the appropriate model to probe next  %</span>
0933     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0934     [hnmin,ii]=max(hn);         <span class="comment">% chose model with the biggest expected decrease</span>
0935     xx=xn(ii);
0936     m=mq;
0937     v=vq;
0938     <span class="keyword">if</span> nr(11)    <span class="comment">% if using log-slope</span>
0939         m(2,:,:)=exp(m(2,:,:));  <span class="comment">% convert to real slope</span>
0940         v(2,:)=v(2,:).*m(2,:,1); <span class="comment">% correct the covariance</span>
0941         v(3,:)=v(3,:).*m(2,:,1).^2; <span class="comment">% and the slope variance</span>
0942     <span class="keyword">end</span>
0943 <span class="keyword">else</span>
0944     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0945     <span class="comment">%     iq=0 means output  model parameters and probe results   %</span>
0946     <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0947     xx=pr;
0948     ii=nr;
0949     m=res(1:nres,:);
0950     <span class="keyword">if</span> ~nargout                     <span class="comment">% print them if no output arguments to call</span>
0951         fid=LOG+(LOG==0);           <span class="comment">% use standard out if LOG==0</span>
0952         pdesc={<span class="string">'Threshold'</span>,<span class="string">'Miss prob'</span>,<span class="string">'Guess Prob'</span>,<span class="string">'Min SNR'</span>,<span class="string">'Max SNR'</span>,<span class="string">'Min Slope'</span>,<span class="string">'Max Slope'</span>};
0953         fprintf(LOG,<span class="string">'\n*********\nModel-specific Parameters\n'</span>);
0954         <span class="keyword">for</span> i=1:7
0955             fprintf(LOG,<span class="string">'%4d) %s: '</span>,i,pdesc{i});
0956             <span class="keyword">if</span> size(pr,2)&gt;1
0957                 fprintf(LOG,<span class="string">'%.5g, '</span>,pr(i,1:size(pr,2)-1));
0958             <span class="keyword">end</span>
0959             fprintf(LOG,<span class="string">'%.5g\n'</span>,pr(i,size(pr,2)));
0960         <span class="keyword">end</span>
0961         qdesc={<span class="string">'nx  SNR values in PDF'</span>, <span class="keyword">...</span>
0962             <span class="string">'ns  Slope values in PDF'</span>, <span class="keyword">...</span>
0963             <span class="string">'nh  Probe SNR values to evaluate'</span>, <span class="keyword">...</span>
0964             <span class="string">'cs  Weighting of slope relative to SRT in cost function'</span>, <span class="keyword">...</span>
0965             <span class="string">'dh  Min step size in dB for probe SNRs'</span>, <span class="keyword">...</span>
0966             <span class="string">'sl  Min slope at threshold'</span>, <span class="keyword">...</span>
0967             <span class="string">'kp  Std deviations of the pdfs to keep'</span>, <span class="keyword">...</span>
0968             <span class="string">'hg  Amount to grow expected gains in ni trials'</span>, <span class="keyword">...</span>
0969             <span class="string">'cf  Cost function'</span>, <span class="keyword">...</span>
0970             <span class="string">'pm  Psychometric model'</span>, <span class="keyword">...</span>
0971             <span class="string">'lg  Use log slope as parameter'</span>, <span class="keyword">...</span>
0972             <span class="string">'pm  Prior std devs in semi-width'</span>, <span class="keyword">...</span>
0973             <span class="string">'pf  Integrated probability floor'</span>, <span class="keyword">...</span>
0974             <span class="string">'ts  Number of std devs to explore'</span>, <span class="keyword">...</span>
0975             <span class="string">'dp  Max shift in probe value'</span>, <span class="keyword">...</span>
0976             <span class="string">'it  Grid interpolation threshold'</span>, <span class="keyword">...</span>
0977             <span class="string">'at  Grid boundary tolerance'</span>, <span class="keyword">...</span>
0978             <span class="string">'la  Look 2-ahead when choosing probe'</span>, <span class="keyword">...</span>
0979             <span class="string">'op  Oulier probability threshold'</span>, <span class="keyword">...</span>
0980             <span class="string">'rx  min grid shrink factor'</span>};
0981         qoptf=[9 10]; <span class="comment">% fields with options</span>
0982         qoptval={<span class="string">'variance'</span>,<span class="string">'entropy'</span>; <span class="keyword">...</span>
0983             <span class="string">'logistic'</span>,<span class="string">'cumulative gaussian'</span>};
0984         fprintf(LOG,<span class="string">'\nShared Parameters\n'</span>);
0985         <span class="keyword">for</span> i=1:length(nr)
0986             fprintf(LOG,<span class="string">'%4d) %s: '</span>,i,qdesc{i});
0987             j=find(qoptf==i,1);
0988             <span class="keyword">if</span> numel(j)
0989                 fprintf(LOG,<span class="string">'%d=%s\n'</span>,nr(i),qoptval{j,nr(i)});
0990             <span class="keyword">else</span>
0991                 fprintf(LOG,<span class="string">'%.5g\n'</span>,nr(i));
0992             <span class="keyword">end</span>
0993         <span class="keyword">end</span>
0994         fprintf(LOG,<span class="string">'\n'</span>);
0995     <span class="keyword">end</span>
0996 <span class="keyword">end</span></pre></div>
<hr><address>Generated by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>